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ABSTRACT 



Context. The Sun is a magnetic star whose magnetism and cyclic activity is linked to the existence of an internal dynamo. 
Aims. We aim to understand the establishment of the solar magnetic 22-yr cycle, its associated butterfly diagram and field parity 
selection through numerical simulations of the solar global dynamo. Inspired by recent observations and 3D simulations that both 
exhibit multicellular flows in the solar convection zone, we seek to characterise the influence of various profiles of circulation on the 
behaviour of solar mean-field dynamo models. We focus our study on a number of specific points: the role played by these flows 
in setting the cycle period and the shape of the butterfly diagram and their influence on the magnetic field parity selection, namely 
the field parity switching from an antisymmetric, dipolar field configuration to a symmetric, mostly quadrupolar one, that has been 
discussed by several authors in the recent literature. 

Methods. We are using 2-D mean field flux transport Babcock-Leighton numerical models in which we test several types of merid- 
ional flows: 1 large single cell, 2 cells in radius and 4 cells per hemisphere. 

Results. We confirm that adding cells in latitude tends to speed up the dynamo cycle whereas adding cells in radius more than triples 
the period. We find that the cycle period in the four cells model is less sensitive to the flow speed than in the other simpler meridional 
circulation profiles studied. Moreover, our studies show that adding cells in radius or in latitude seems to favour the parity switching 
to a quadrupolar solution. 

Conclusions. According to our numerical models, the observed 22-yr cycle and dipolar parity is easily reproduced by models includ- 
ing multicellular meridional flows. On the contrary, the resulting butterfly diagram and phase relationship between the toroidal and 
poloidal fields are aff'ected to a point where it is unlikely that such multicellular meridional flows persist for a long period of time 
inside the Sun, without having to reconsider the model itself. 

Key words. Sun: magnetic fields - Sun: activity - Sun: interior - Methods: numerical 



1. Introduction 

The Sun possesses striking magnetic and dynamical properties, 
such as its turbulent convective envelope, large-scale surface 
diflTerential rotation, 22-yr cycle of magnetic activity, butterfly 
diagram of sunspot emergence, hot corona, etc. (Stix 120021) . 
Understanding how the physical processes operating in the so- 
lar turbulent plasma nonlinearly interact to yield this wide range 
of dynamical phenomena is very challenging. One successful 
and powerful approach is to rely on multi-dimensional magne- 
tohydrodynamics (MHD) numerical simulations. Today, despite 
tremendous advances in building powerful supercomputers, it is 
still not possible to compute a fully integrated 3-D MHD model 
of the Sun starting from its core up to its corona. One is thus 
forced to study individually complementary pieces of the full 
solar MHD puzzle and to progressively incorporate them in a 
more nonlinearly coupled model. One important characteristic 
of the Sun that needs to be understood is the origin of its mag- 
netic activity because it has direct societal impact by impairing 
satellites, damaging electric power grids, interfering with high 
frequency radio communications and radars. It is currently be- 
lieved that the solar magnetism is linked to an internal dynamo 
(Parker 1955 ). More precisely the Sun is the seat of both a small 
scale and irregular dynamo and a large scale and cyclic dynamo 
that generate and maintain its magnetic field and lead to the var- 
ious magnetic phenomena observed at its surface (Parker ] 19931 
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Cattaneo & Hughes 12001 [ Ossendrijver 120031) . Developing nu- 
merical models of the solar dynamo has thus been a very active 
field of research. This has mainly involved two types of numeri- 
cal experiments: 

- kinematic solar dynamo models that solve only the induc- 
tion equation in its mean field approximation and assume the 
velocity field as given (Steenbeck & Krause |1969t Roberts 
fT972l Stix[T97g: MoflTat 1978; Krause & Ra dler [19801 see 
Charbonneau [2005 and Solanki et al. 120061 for recent re- 
views). These models rely on the parametrization of two im- 
portant eflfects that are thought to be at the origin of the solar 
global dynamo, the a and Q effects. They provide a useful 
and fast tool to model the solar 22-yr magnetic cycle and 
its associated butterfly diagram since no feedback from the 
Laplace force on the motion is accounted for. 

- or dynamical solar dynamo models that solve explicitly 
the fu ll set of MHD equations (Oilma n [T9831 Glatzmaier 
[T9851 Cattaneo 1999; Brun et aL l2QQ4b . These models self- 
consistently compute all the physical processes in three di- 
mensions allowing significant progress to be made on the 
intricate interactions operating in a turbulent magnetized 
plasma. The cost of 3D models and the large number of de- 
grees of freedom needed to model the whole Sun make it 
difl&cult, as of today, to provide quantitative predictions such 
as the cycle period. 

Clearly, both approaches are complementary and are needed to 
better understand the magnetic solar activity. Since the original 
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ideas of Parker regarding the operation of a hydromagnetic dy- 
namo in the Sun, many articles have been written to improve our 
understanding of this subtle physical process. In the late 70's, 
solar dynamo models were relying on a cylindrical differential 
rotation profile and an of-eff'ect linked to non-reflexion symmet- 
ric motions within the turbulent and rotating solar convection 
zone (CZ), the so called a - dynamo. In such dynamo mod- 
els, the product of a and dfl/dr must be negative in the northern 
hemisphere in order to obtain an equatorward butterfly diagram 
(Yoshimura 1975). However, these distributed a - dynamo 
models have since been discarded for two main reasons: first 
the inversion in the mid 80 's of the internal solar rotation pro- 
file (Brown et al. 1989; Thompson et al. 2003 ) showed a conical 
diflTerential rotation profile in the convection zone (dQ/dr ^ 0) 
rather than a cylindrical profile. Secondly, it was demonstrated 
that strong magnetic fields could significantly reduce the ef- 
ficiency of the Qf-eflTect, in a phenomenon called a-quenching 
(Ossendrijver 2003). In a landmark paper, Parker (1993) pro- 
posed the segregation of sites of generation of the poloidal field 
on the one hands with that of the toroidal field in the other hand, 
in what is now called the interface dynamo. He was encouraged 
by the latest helioseismic inversions which indicated the exis- 
tence of a swift transition from the diflTerential rotation of the 
solar convection zone to an inner solid body rotation in the ra- 
diative interior, i.e the tachocline (Spiegel & Zahn 1992). In the 
late 90's, Charbonneau & Mc Gregor (1997), were the first to 
incorporate all the ingredients of the modern interface dynamo: 
a solar-like (conical) diflTerential rotation -h a tachocline, a sepa- 
rate site of generation of poloidal field (in the convection zone) 
vs the toroidal field (in the tachocline). They showed that with 
this new solar dynamo model, the 22-yr cycle period, the but- 
terfly diagram, the phase relationship between the poloidal and 
toroidal fields and the field parity can be reproduced. However, 
these models do not include meridional circulation (MC). To 
address this issue, Dikpati & Charbonneau (I1999D computed 
Babcock-Leight on (BL ) models (Babcock 1962 Leighton 1969; 
Choudhuri et al. 11995b with a solar-like Q profile and an unicel- 
lular meridional flow. They showed that a solar dynamo model 
based on this so-called Babcock-Leighton flux transport dynamo 
could also be successful at reproducing most of the solar global 
magnetic properties. In this model, the meridional circulation 
transports the poloidal field from the surface, where it appears 
through the twisted nature of the solar active regions, to the bot- 
tom of the convection zone where it is transformed into a toroidal 
field in the tachocline. This meridional circulation thus plays a 
major role in the behaviour of BL flux transport dynamo models. 
It is then important to understand its origin and structure in the 
Sun. 

An analysis of the governing equations tells us that mean 
meridional flows arise from a combination of buoyancy forces, 
Reynolds stresses, latitudinal pressure gradients and Coriolis 
forces acting on the mean zonal flow (diflTerential rotation) 
(Miesch 120051) . The competition of these physical processes 
make it diflScult to anticipate the meridional flow profile. Inside 
the solar envelope, this flow is much weaker than the diflTerential 
rotation, making it relatively difl&cult to measure. Furthermore, 
although it can in principle be probed by global helioseismology, 
its eflTect on global acoustic waves is weak and difliicult to distin- 
guish from rotational and magnetic eflTects. Thus, we must cur- 
rently rely on surface measurements and local helioseismology. 
The 15 m.s"^ poleward flow observed at the surface (Hathaway 
119961) has been confirmed by local helioseismology with great 
accuracy down to r/R^ = 0.95 (Haber et al. 120021) and some 
attempts have been made to probe the MC down to r/RQ = 0.85 



(Giles et al. [T997l Schou & Bogart [T998l Braun & Fan [T998l) but 
the pattern and localisation of the equatorward return flow is still 
not well established. Today, the favoured solar dynamo models 
are of flux transport type, assuming both a source of poloidal 
field at the surface (a BL source term) and at the bottom (a- 
eflTect like) (Bonanno et al.[2QQl Dikpati et al. [20041 Kuker et al. 
2001 ; Chatterjee et al. 2004). In particular, recently these models 
have been successful at reproducing a series of solar cycle and 
even for starting to predict the next/starting solar cycle (cycle 
24) (Dikpati & GilmangOOS). 

In this paper we will follow the kinematic approach, by com- 
puting 2-D axisymmetric mean field solar dynamo models of the 
flux transport BL type. We seek to answer the simple follow- 
ing questions: What is the role of meridional flows in setting 
the solar cycle period and butterfly diagram? Can the presence 
of multicellular meridional flows lead to variations of the gen- 
eral properties of the solar activity? The motivation behind these 
questions is that both observational evidence via local helioseis- 
mology technics (Haber et al. 120021) and 3-D MHD numerical 
models as described above (Miesch et al. '2000; Brun & Toomre 
2002; Brun et al. 2004; Browning et al. 2006 ) exhibit multicellu- 
lar flow both in radius and latitude. If such permanent multicel- 
lular flow were indeed acting continuously in the Sun, it is likely 
that it will lead to a diflferent solar global dynamo model since to- 
day most models rely on a single monolithic meridional flow to 
transport poloidal field from the surface down to the tachocline 
at the base of the solar convection zone. 

The paper is organized as follows: in Sect. 2, we present the 
mean field induction equation and the ingredients of the model, 
in Sect. 3, we discuss the results of our study, mainly the ef- 
fect of introducing many meridional cells both in latitude and 
radius in the model. In Sect. 4 we discuss the influence of the 
more complex meridional flow in setting the field parity (i.e ei- 
ther dipolar or quadrupolar) and we conclude in Sect. 5. Finally, 
the numerical techniques used to solve the induction equation 
and the boundary conditions introduced to compute the tempo- 
ral evolution of our solar dynamo models are presented in the 
appendix. 

2. Setting the solar dynamo model 

2.1. Mean field equations 

To model the solar dynamo, we use the hydromagnetic induction 
equation, governing the evolution of the large scale magnetic 
field B in response to advection by a flow field v and resistive 
dissipation. 

— = V X (v X B) - V X (77V X B) 
ot 

where 77 is the eflfective magnetic diflTusivity. 

Working in spherical coordinates and under the assumption 
of axisymmetry, we write the total mean magnetic field B and 
the velocity field v as: 

B(r, 6/, = V X (A^(r, 6, t)e^) B^(r, 6, t)e^ 

v(r, 0) = Vp(r, 0)-\-r sin Q(r, 0)e^ 

Note that our velocity field is time-independant since we 
will not assume any fluctuations in time of the diflferential ro- 
tation Q or of the meridional circulation Vp. Reintroducing this 
poloidal/toroidal decomposition of the field in the mean induc- 
tion equation, we get two coupled partial diflTerential equations. 
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one involving the poloidal potential and the other concerning 
the toroidal field 



1 



(1) 



dt 



1 1 d(TnB,) d{nln,) B, 



TU 



dr 



TU 



(2) 



where = r sin 6, % is the turbulent magnetic diff'usivity 
(diff'usivity in the convective zone), Vp the flow in the meridional 
plane (i.e. the meridional circulation), Q the difl'erential rotation. 
The break of axisymmetry needed to circumvent Cowling's anti- 
dynamo theorem comes from the addition of a term S (r, 6, 5^) 
in Eq.©, representing the BL surface source term for poloidal 
field. In order to write these equations in a dimensionless form, 
we choose as length scale the solar radius Rq and as time 
scale the diff'usion time Rq/tji based on the envelope diff'usiv- 
ity T/f This leads to the appearance of three control parame- 
ters Cq = QoR^/rju Cs = soRq/tji and Rq = vo^o/^t where 
^o^sq, vo are respectively the amplitude of the diff'erential rota- 
tion, of the surface source term and of the meridional flow. 

Equations [T] and [2] are solved in an annular meridional cut 
with the colatitude 6 e [0,7r] and the radius r g [0.6, 1]7?q 
i.e from slightly below the tachocline (r = O.IRq) up to the 
solar surface, using a finite element method (STELEM code) 
which was validated thanks to an international dynamo bench- 
mark (Jouve et al. .2007i see the appendix for more details on the 
numerical technique). At ^ = and = n boundaries, both A^ 
and B(p are set to and at r = 0.67?q, both A^ and B(p are set to 
. At the upper boundary, we smoothly match our solution to 
an external potential field, i.e. we have vacuum for r > Rq. As 
initial conditions we are setting a confined dipolar field config- 
uration, i.e the poloidal field is set to sinO/r^ in the convective 
zone and to below the tachocline whereas the toroidal field is 
set to everywhere. 



2.2. The physical ingredients 

The model "ingredients" are basically those used by Dikpati & 
Charbonneau (1999). The rotation profile is a representation of 
that deduced from helioseismic inversions, assuming a solid ro- 
tation below 0.657?© and a diff'erential rotation above the in- 
terface. With this profile, the radial shear is maximal in the 
tachocline: 



Q(r, ^) = Qc + ^ (^Eq + ci2 cos^ 6-\-a4 cos^ 6 - Qc) 



1 -herf 2 



with I^Eq = 1, = 0.93944, = 0.1 Rq, di = 0.05, 
a2 = -0.136076 and = -0.145713. 

In BL flux transport dynamo models, the poloidal field owes 
its origin to the twist of the magnetic field emerging at the solar 
surface. Thus, the source has to be confined in a thin layer just 
below the surface and as the process is fundamentally non-local, 
the source term depends on the variation of B(p at the base of the 



convection zone. Moreover, a quenching term is introduced to 
prevent the magnetic energy from growing exponentially: 



S(r,e,B^) 



r — rn 
l+erf(— ^) 



1 + 



l-erf(^) 

d2 



cos0sm6B(p(rc, 0, t) 



where r2 = 0.95, = 0.01, = 10^ 

We assume that the net diff'usivity in the envelope r] is dom- 
inated by its turbulent contribution whereas in the stable zone, 
the value of the diff'usivity has to be much weaker (we have 

« Vi)- We smoothly match the two diff'erent constant val- 
ues thanks to an error function which enables us to quickly and 
continuously transit from 7/^ = 10^ cm^.s"^ to 7/t which is a vari- 
able parameter in our computations. It gives us the diflfusivity 
function below: 



1 -herf 2 



1 = ^ + 1 

We have now set up a detailed model of the global solar dy- 
namo, using the framework of mean field theory. One of the key 
ingredients of this kind of models is the meridional circulation. 
We now investigate the influence of complex flows on the solar 
dynamo and its global properties. 



3. Influence of meridional circulation on magnetic 
cycles 

3.1. Our reference unicellular model 

We first compute a model where we assume one large single 
meridional cell per hemisphere which we will consider as the 
reference model. The components of the meridional circulation 
are those used in Van Ballegooijen & Choudhuri (1988) which 
defines a steady circulation pattern, symmetric with respect to 
the equator, with a single flow cell per hemisphere directed pole- 
ward at the surface and allowed to penetrate a little below the 
base of the CZ, where it is equatorward. With this typical model, 
we are able to reproduce several aspects of the solar cycle, no- 
tably its period of approximately 20 years, a strong equatorward 
branch for toroidal field restricted to low latitudes, a phase shift 
of tt/ 2 between the surface polar field and the deep toroidal field, 
so that the polar field changes its polarity from negative to pos- 
itive when the toroidal field is positive and maximal in inten- 
sity near the equator (Fig.[T]). Moreover, the strong equatorward 
branch for the toroidal field is the signature of the drag of the 
toroidal field by equatorward MC at the base of the convection 
zone and thus clearly shows the dominating eff'ect of field ad- 
vection over diff'usion. Indeed, a least square fit indicates that 
the cycle period T strongly depends on the meridional flow am- 
plitude, i.e. T oc v~^-^^. 

One could wonder why we are seeking to improve and mod- 
ify the reference model given its relatively good agreements with 
observations. There are in fact several reasons. 

First, a 26-year interval studied by Snodgrass and Dailey 
( 1996) exhibited large temporal variations in the meridional flow 
amplitude. Indeed they found that even if the latitudinal flow 
peaked at about 15m.s"^ in average, MC could achieve ampli- 
tudes as large as 50m.s"^. Given the strong dependance of the 
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Fig. 1. Reference case: butterfly diagram (time-latitude cut at 
r = est) of the unicellular model with vq = 643cm.s"\ = 
20cm.s"^ and 7/t = 5.10^^ cm^.s"^. The contours of 5^ (upper 
panel) are plotted at the base of the convection zone and By 
(lower panel) is taken at the surface. Contours are logarithmi- 
cally spaced with 2 contours covering a decade in field strength 
and red colours represent positive values of the field. The vertical 
dashed line corresponds to the epoch of reversal of toroidal field, 
the plain line correspond to the epoch of reversal of poloidal field 
at the poles from negative to positive polarity and the dash-dotted 
line corresponds to the positive maximum of toroidal field near 
the equator. 



cycle period on the amplitude of the flow (if we triple the ve- 
locity amplitude we reduce the period by about one-half), we 
can thus wonder if a systematic period of about 22 years can be 
conserved in this context of temporally varying flows. A study 
of the impact of stochasticity in such BL models (Charbonneau 
& Dikpati 2000) yet showed that these solar cycle models were 
quite robust to stochastic variations of the meridional circulation. 

Another argument against this type of BL models comes 
from Dikpati & Charbonneau (1999) who showed that even if 
the configuration of the toroidal field seems to fit the obser- 
vations quite well particularly concerning the strong equator- 
ward branch, a relatively strong toroidal field (10^ G) is also 
present at all latitudes. This strong field existing at all latitudes 
could be significantly decreased by imposing a lower thresh- 
old for quenching in the surface source term that would prevent 
toroidal flux tubes that are too weak in intensity to rise through 
the CZ and thus participate in the regeneration of the poloidal 
field (Charbonneau et al. 2005 ). 

Finally, Dikpati & Oilman raised in 2001 a major concern 
about the BL flux transport model, concerning the symmetry of 
the magnetic field with respect to the equator. They claim that, in 
the particular range of parameters they are using to get a solar- 
like period, the pure BL flux transport model fails to reproduce 
the persistent antisymmetry of the toroidal field and that what- 
ever the magnetic initial conditions imposed, this model would 
always end up giving a quadrupolar configuration, which we do 
not currently observe in the Sun (see Sect.|4l). 

Thus the single cell pure BL model does not seem fully satis- 
factory and needs to be improved. Moreover, both observations 



by Haber et al. (I2002D and 3D simulations by Brun et al. (I2004D 
show multiple cells circulation in the CZ and modulation of the 
MC with magnetic fields. Dikpati et al. (2004) and Bonanno et 
al. (2005 ) computed dynamo models including 2 cells in lati- 
tude per hemisphere. Dikpati, with a BL source term, found that 
these additional cells tended to decrease the cycle period and 
Bonanno, with a distibuted CK-eflTect model, found that the global 
pattern of the meridional circulation could strongly influence 
the location of the dynamo action in the advection-dominated 
regime. Consequently, we would like to verify the influence of 
even more complex multicellular (both in radius and in latitude) 
meridional flow on the cycle, on the butterfly diagram, on the 
phase relationship between the poloidal and the toroidal parts of 
the magnetic field. 

3.2. Solar dynamo models with additional cells in the 
mendional circulation 

We focus here on two cases, the case with 2 cells in radius and 
the case with 2 cells in radius and 2 in latitude. We will not deal 
with the 2 latitudinal cells model, as this was already treated 
by Dikpati et al. (12004 ) and Bonanno et al. (12005b . The main 
characteristics of the diflTerent cases studied are summarized in 
Tablell] 

To get a multicellular flow, we write the stream function as 
a product of Chebyshev polynomials in radius and of Legendre 
polynomials in latitude. Through the following equations: pvr = 
p- ^ and pve = - ^ ^ , with z = r, x = - cos 6 and p = 

1/z, which ensure that V.(pv) = 0, we easily deduce the shape 
of the meridional flow components from the polynomial stream 
function. 

3.2.1 . Case 1 : 2 cells in radius per hemisphere 

For this case, all the ingredients are kept identical to the ref- 
erence model but the meridional flow is modified. We set the 
stream function i// to: 

jA(x, z) = Kiziz - 0.6)(500(z - 0.8)^ - 20(z - 0.8))(x^ - x) 

Ki is a normalization factor, i.e. it is chosen so that V6i/vo = 1 
at the solar surface and at a latitude of 45° (see Fig.O. 

In Fig. [3l we represent the butterfly diagram of case 1 with 
the parameters used in the reference unicellular model. For this 
model the cycle lasts 84.6 years, more than 3 times longer. The 
increase of the period comes from the fact that the magnetic flux 
is not transported from the surface to the interface as fast as it 
was in the unicellular model. This is a direct consequence of 
the presence of a return flow at mid depth. The two source re- 
gions (the surface for the poloidal field and the tachocline for 
the toroidal field) are thus not linked as directly as they were in 
the reference model. This leads to a slower regeneration of the 
toroidal field from poloidal field and vice versa. 

We also see that since the meridional flow is directed pole- 
ward at the base of the convective zone, we get a very strong 
poleward branch for the toroidal field. Moreover, the last panel 
of Fig. [2] shows that the latitudinal velocity intensity is 3 times 
higher at mid-depth (where it is equatorward) than at the base of 
the CZ (where it is poleward). As a consequence, the field is ad- 
vected 3 times faster to the equator at mid-depth than to the poles 
at the base of the convection zone and this explains the strong 
domination over time of the poleward against the equatorward 
branch seen on the upper panel of Fig.[3l This figure also shows 
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Table 1. Summary of the 5 different cases and associated parameters. The last column indicates the period in years for each case. 





Resolution 


Time step 


vo 




So 


Cycle Period 




X 




(cm.s~^) 


(cm^.s~^) 


(cm.s~^) 


(Yrs) 


Reference case 


128^ 


7.2 10"' 


643 


5.10^^ 


20 


21.8 


Case la 


256 X 128 


4.53 10"' 


643 


5.10^^ 


20 


84.6 


Case lb 


256 X 128 


4.53 10-^ 


1916 


1.4910^^ 


20 


22.4 


Case 2a 


256 X 128 


4.53 10-^ 


643 


5.10^^ 


20 


44.7 


Case 2b 


128^ 


7.2 10-^ 


1071 


1.510^^ 


20 


22.4 



2 radial cells 
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Fig. 2. Stream function and components of the meridional flow 
multiplied by vq = 643 cm.s"^ for the 2 radial cells model. 



the existence of an equatorward branch between and 30° as re- 
quested by observations. A small amount of poloidal field is thus 
advected towards the equator even though the flow is poleward 
in this region. This is a direct consequence of the non-locality of 
our surface source-term. Indeed, at these latitudes and between 
0.737?© and 0.947?©, the toroidal field is advected toward the 
equator by the MC flow. Both through advection and difl'usion, 
this equatorward-migrating toroidal structure is transported in- 
ward near 0.77?© and as the poloidal source term is linked to the 
toroidal field at this radius, the poloidal field ends up drifting 
equatorward. This figure also shows the appearance of smaller 
scale structures in the radial field at the surface. In particular, 
we note a small equatorward branch (about 80 times less intense 
than the value near the pole) in a very narrow band around the 
equator (between -20° and 20° in latitude) which is of opposite 
polarity than that of the present cycle. This branch is the rem- 
nant of a small amount of field of the preceeding cycle which 
was driven back up to the surface by the upper cell which creates 
an upflow at mid-depth near the equator (see the radial velocity 
profile on Fig. [2]). 





200 
Time (Years) 



Fig. 3. Case la: butterfly diagram (time-latitude cut at r = est) 
of case 1 with vq = 643 cm.s"^. The format is the same as Fig.[T] 



The dependance of the period of this model on variation of 
the magnetic Reynolds number and thus of the velocity am- 
plitude is very strong in comparison to the unicellular model. 
In this case, we have the following dependance for the period: 
T oc v~^-^^. The strong dependance of the cycle period on the 
MC amplitude suggests that it would be easy to recover a solar 
period of about 22 years, only by increasing the amplitude of 
the meridional flow, keeping the other parameters constant. So 
the maximum latitudinal velocity vq needed to get a 22-yr cycle 
period keeping = 20cm.s"\77t = 5 x lO^^cm^.s"^ would be 
about 2500 cm.s"^. However, if we only increase the MC ampli- 
tude, we lose the antisymmetry of the toroidal field with respect 
to the equator observed in the Sun (see Sect. |4]). Thus, to keep 
the correct dipolar parity for this 2 radial cell model, we need 
to increase both the MC amplitude and the magnetic diff'usivity, 
hence, we get a solar-like parity 22-yr model with the follow- 
ing parameters: ^0 = 20cm.s"\7/t = 1.49 X lO^^cm^.s"^ and 
1916 cm. s"^ In Fig.Hl we represent the butterfly diagram for this 
22-yr cycle case which also exhibits the phase relationship be- 
tween the poloidal and toroidal fields. We maintain the equator- 
ward branch for the toroidal field between and 30°. Moreover, 
the radial field evolution is smoothened by the increased diff'u- 
sivity and we thus see much fewer small structures in the lower 
panel of Fig. (H We nevertheless note that we keep very strong 
values for the polar field at the surface, which was also the case 
for the reference model. 

Adding a new cell in radius modifies the magnetic advective 
path and thus the link between the two source regions (the sur- 
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Time 



Fig. 4. Case lb: butterfly diagram (time-latitude cut at r = est) 
of case lb with vq = 1916 cm.s"^ and 7]t = 1.49 x 10^^ cm^.s"^. 
The format is the same as Fig. [H 



face and the base of the convection zone). As a consequence, 
the time-delay between the reversal of the polar field at the sur- 
face and the maximum of toroidal field at the base of the CZ is 
also modified. Indeed, the phase shift between the 2 components 
of the magnetic field is here about 7r/3: we observe that as the 
poloidal field reverses at the pole, the toroidal field has not yet 
reached its maximum, thus lagging the poloidal field. 

In Fig.[5l we show the field evolution in the meridional plane 
of the 22-year cycle model. We see that the field configuration 
tends to follow the complex nature of the MC. Indeed, we clearly 
see that we are in the advection-dominated regime as the major 
field concentration areas follow the meridional flow streamlines. 
Figure [2l indicates that we have an upper cell with a poleward 
flow very concentrated in a thin layer near the surface which 
we recover on the magnetic patterns especially on panels c) and 
d) where a small part of the toroidal field is being advected to- 
wards the pole near the surface. Between 0.73 and 0.94Rq (i.e in 
more than 60% of the CZ), where the flow is equatorward, most 
toroidal field of a chosen polarity is advected toward the equator 
and amplified by the latitudinal shear of the poloidal field. As it 
reaches the equator, it splits in two parts as we see on panel J), 
one being redirected towards the surface where it will be driven 
in the direction of the pole by the top of the upper meridional 
cell and the other part, containing most of the toroidal field, be- 
ing advected towards the base of the CZ. As the toroidal field 
reaches the base of the CZ, the poleward flow advects the field 
towards the pole [panels e) and /)] where it is amplified by the 
radial shear of poloidal field which at the same time makes the 
opposite polarity of the preceeding cycle decay away. We can 
also see that the poleward branch of the toroidal field of one cy- 
cle (cycle n) is "pushed" towards the pole both by the field of 
cycle n - I and cycle ^ - 2, all present at the base of the CZ at 
the same period of time. Panels J), e) and /) show that the fields 
of cycle n and ^ - 2, of the same polarity, even reconnect with 
each other at a given time and stay connected during almost half 
a magnetic cycle, before being split back by the diverging cells 
of meridional flow. 




Fig. 5. Case lb: temporal evolution of the poloidal potential (left 
panel) and the toroidal field contours (right panel) in a merid- 
ional plane for case 1 during half a magnetic cycle. The blue 
contours indicate an anticlockwise orientation for Bp and a nega- 
tive orientation for 5^ (i.e the field is directed towards the reader 
in this case) and the red ones clockwise orientation for Bp and a 
positive orientation for B^p (directed away from the reader). 
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We thus see that having many cells in radius impacts signif- 
icantly BL models. We now turn to studying the coupled effect 
of having two cells in radius and latitude. 

3.2.2. Case 2: 2 cells in radius and 2 cells in latitude per 
hemisphere 

In this model, the stream function </r is a product of polynomials 
of higher order (degree 5 in z and x) that we have to multiply by 
a function of x which enables us to choose the relative velocity 
amplitude in each cell. We set the 2 top cells to the same maxi- 
mum value which implies that i// has the following expression: 



i//(x,z) = i^2z(z-0.65)(500(z -0.825)^- 500 X (0.175)^ 
X (z- 0.825)) X (Ix^ - lOx^ -h 3x) x (1 - x^f^^-^^^ 

if z> 0.65 and otherwise (see fig[6|. 
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Fig. 6. Stream function and components of the meridional flow 
multiplied by vq = 643 cm.s"^ for case 2. 



For the 4 cells model, the butterfly diagram corresponding to 
the exact same parameters rjt, sq and vq as the unicellular model 
is shown in Fig. [71 The most obvious property of this model is 
that it can sustain two magnetic cycles, one near the equator and 
the other at high latitudes with significantly diff'erent periods. In 
both equator and polar branches, the cycle period is strongly in- 
creased, up to 44.7 years near the equator and reaching a period 
of 124 years at high latitudes. This behaviour is due to the pres- 
ence of 2 cells in radius which significantly increases the period 
probably because, as in case 1, the magnetic flux is not trans- 
ported from the surface to the interface as fast as it was in the 
unicellular model because a return flow is present at mid depth. 
However, adding cells in latitude decreases the time for the fluid 



Fig. 7. Case 2a: butterfly diagram (time-latitude cut at r = est) 
of case 2 with vq = 643 cm.s"^. The format is the same as Fig.[T] 
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Fig. 8. Case 2b: Butterfly diagram (time-latitude cut at r = est) 
and field phase relation of case 2 with a period close to the solar 
period. 



to travel along the "conveyor belt" for the low-latitude region, 
hence the flux is transported faster from the surface to the base 
of the convection zone and thus the regeneration of each com- 
ponent of the magnetic field is faster, as shown by Dikpati et al. 
( 2004). Here, by having cells both in radius and in latitude, we 
get a cycle faster than the two radial cells case but still slower 
than the unicellular case. Consequently, the influence of hav- 
ing several radial cells seems to be much stronger than that of 
adding cells in latitude. We can see on this butterfly diagram 
that the patterns are quite complex in the time-latitude plane. We 
clearly see on the toroidal field at the base of the CZ and on the 
poloidal field at the surface the imprint of the 2 counter-cells at 
45°. Once again the magnetic field behaviour strongly depends 
on the direction of the flow. Indeed, at the base of the CZ, the 
toroidal field at high latitudes (where the flow is equatorward) is 
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drifting from 60° to 45° where it encounters the counter latitudi- 
nal cell. At low latitudes (where the flow is poleward), the field 
is advected polewards from the equator to the zone of vanish- 
ing V6I (i.e. around 45°). For the radial field, the patterns are also 
very intricate. We again find the small equatorward branch as in 
case 1 . We also clearly see the generation of new poloidal field 
structures just above 50° embedded in the structure of opposite 
polarity. It is then moving down to 45°, meeting a branch of the 
same polarity coming from the preceeding cycle and from the 
equator. The main part of the radial field at the surface is then 
advected towards the pole, creating a significant polar branch in 
the high-latitudes part. One interesting feature of this model is 
that the intensity of the polar field is diminished by a factor 10 
compared to the unicell model, which could be in better agree- 
ment with the solar observations. 

We tested the sensitivity of the multicellular model's charac- 
teristics to variations of the physical parameters. Using a least 
square fit to get the exponents of each parameter, we note that 
the dependance of the cycle period at low latitudes on ^o, vq and 
rii is as follows: 

The dependance of the period of this model to a variation of 
the velocity amplitude is reduced in comparison to the previous 
models. The cycle period will thus be less disturbed by tempo- 
ral fluctuations of the MC amplitude, which is a very attractive 
feature of this model. For any value of ^o, we expect that an in- 
crease in the intensity of the meridional flow enables the field to 
travel faster along the 'conveyor belt' so that both poloidal and 
toroidal fields head faster towards their reversals. This is why we 
have a negative dependance of the period on vq. 

As soon as the Reynolds number becomes too high (Rq above 
800), i.e. when the strength of the meridional flow is increased, 
77t remaining constant, a strong polar branch with a longer period 
is appearing. This property is due to the advection by meridional 
flow of the magnetic field in the whole convection zone dom- 
inated by its poleward component (see below). When the field 
reaches the base of the CZ, the strong toroidal structure is thus 
concentrated near the pole, trapped in the slowly moving merid- 
ional cell at high latitudes. As we do not get this feature of two 
coexisting branches of different periodicity in the Sun, we seek 
to recover a unique cycle period, taking into account the depen- 
dance on parameters. We have seen that it is not sufl&cient to act 
on the strength of the meridional circulation to recover the 22-yr 
cycle period, but as the least square fit shows, turbulent diflTusiv- 
ity plays a crucial role in this multicellular model. 

The strong negative dependance of the cycle period on the 
turbulent magnetic diflTusivity is characteristic of the fact that 
we have multiple cells of meridional flow in each hemisphere. 
The magnetic field follows the configuration of the meridional 
flow as long as it is advected by the circulation but the magnetic 
diflTusivity enables the field to cross the strong velocity gradi- 
ents present at the borders of each cell. Magnetic diflfusivity thus 
provides the field a way to short-circuit the complex "advection 
path" of this model and to allow for a faster link between one 
meridional cell and another. As a consequence, if 7/t is too low, 
the poloidal field created at the surface will cross the strong ve- 
locity gradients that we have at the borders between the vari- 
ous cells in a much longer time and the classical mechanism 
of regeneration of toroidal field from this existing poloidal field 



will be less efl&cient. Consequently, to obtain a 22-yr cycle pe- 
riod for this model, we have to increase the magnetic diflTusiv- 
ity (without of course losing the advection-dominated regime) 
to 1.5 X 10^^ cm^.s"^ (so that the magnetic field is enabled to 
get from one circulation cell to another quickly enough) and in 
agreement with the least- square fit which indicates a neglige- 
able dependance of the period on and a negative dependance 
on the velocity amplitude, we keep the source term intensity to 
^0 = 20 cm.s"\ and slightly increase vq to 1071 cm.s"^. Figured 
shows the butterfly diagram and the phase relationship between 
the poloidal and toroidal fields for this 22-yr period case. The 
first positive result of this model is that no second cycle period 
appears, we get a unique cycle due to the dynamo action in the 
whole CZ. Here again, as the link between the 2 source regions 
of poloidal and toroidal fields is complex, the field relation is 
not that which we observe in the Sun and it is even more per- 
turbed than in case 1. Indeed, the poloidal field here reverses 
approximately where the toroidal field of the opposite polarity 
reverses, meaning that we have a phase shift of about n between 
the 2 components of the field, with the toroidal field leading 
the poloidal field. We note that in this case, unlike case 1, the 
shape of the butterfly diagram is slightly diflTerent from that using 
the parameters of the single cell model. Indeed, increasing the 
magnetic diflTusivity caused the sharp magnetic structures to be 
smoothed, especially for the radial field in which the small equa- 
torward branch of opposite polarity is not visible anymore and 
where the branch at 50° drifting to 45° is thicker and smoother. 

Figure [9] shows the temporal evolution of the magnetic field 
in the meridional plane for the same case 2b. The first appear- 
ance of a new poloidal field occurs at a latitude of about 30°, 
created by the BL source term. The creation of poloidal field at 
such latitudes is a direct consequence of the non-locality of our 
source term whose latitudinal dependance is linked to the latitu- 
dinal dependance of the toroidal field existing at the base of the 
CZ at the same time. This poloidal field structure is then both 
slowly dragged towards the 45° latitude by the MC and ampli- 
fied by the source term which is still active in this region. The 
poloidal field is here latitudinally sheared to create and amplify 
toroidal field [see panels a) and b)] . After the reconnection, MC 
has dragged the field deeper down to regions of both poleward 
and equatorward flow so the field is located at the particular con- 
vergence point between the 4 cells. This inward advection due 
to the meridional flow configuration at 45° is also acting on the 
toroidal field to drive it down to the middle of the CZ and should 
not be mistaken with diflfusive eflTect. It is really the strongly neg- 
ative radial velocity at this latitude (see panel 2 of Fig.© which 
advects the field inward to the point of convergence of the 4 cells. 
We see on panels c), d) and e) that even if some of the field is 
advected towards the equator, the north branch is dominant and 
thus most of the field goes up to the pole, while the older field of 
opposite polarity is being advected to the base of the CZ. This 
poloidal field of opposite polarity (the negative field on panels 
d), e) and /)) is then heading back to the 45° region where it 
meets another poloidal structure of the same polarity. If we look 
at the toroidal field, we see it is created and amplified in the 
whole convection by the latitudinal shear of the diflTerential ro- 
tation. At the same time, 5^ is being advected by the flow, first 
inward to the middle of the CZ and then mainly to the polar re- 
gions. It is finally dragged down to the base of the CZ where 
the field of the preceeding cycle is cancelled by the creation of a 
new field of opposite polarity. Unlike case 1 , the suppression of 
the field of the preceeding cycle seems thus to be mainly due to 
the advection of the present 5^ (which was created before by the 
shear in latitude of the poloidal field) than to the radial shear of 
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Fig. 9. Case 2b: temporal evolution of the poloidal field lines 
with its potential extrapolation (left panel) and the toroidal field 
contours (left panel) in a meridional plane for case 2 for half a 
magnetic cycle. The format is identical to Fig.[5l 



Bp in the tachocline. In this case, only cycle n and n - I seem to 
really interact since the toroidal field of cycle n-2 has already 
completely vanished at the bottom of the CZ when B^p of cycle n 
is being created in the upper part of the CZ. 

This 4 cells model is thus very intriguing for many impor- 
tant solar dynamo properties. The butterfly diagram as well as 
the field lines evolution during a cycle become very complex, 
the field phase relationship is not corresponding anymore to the 
solar observations but the dependance on the MC amplitude is 
reduced and the strength of the polar field is decreased, which 
constitute attractive characteristics of the model. 



4. Parity selection in multiple meridional cells 
dynamo models 

As we said before, Dikpati & Oilman in 2001 showed that with 
a set of parameters they found appropriate to give a solar-like 
solution, their pure BL flux transport model had difficulties re- 
producing the persistent antisymmetry of the toroidal field we 
observe in the Sun. 

Several solutions were proposed to solve this problem, 
Dikpati & Oilman 6200 1) as well as Bonanno et al. (12002b man- 
aged to get rid of this field parity drift by imposing an a-eff'ect 
at the base of the CZ, thus imposing two spatially separated 
source terms for the poloidal field. Another solution was pro- 
posed by Chatterjee et al. (12004b : they keep the regular sur- 
face source term of BL type but they impose a small diff'usivity 
(2.2 X 10^ cm^.s"^) in the overshoot layer to prevent the toroidal 
field from diff'using across the equator and a very large diff'usiv- 
ity (2.4 X lO^^cm^.s"^) for the poloidal field in the convective 
zone to allow diff'usive coupling of the poloidal field between 
the two hemispheres. 

We now seek to characterise the influence of multicellular 
flows on parity selection and we consider its sensitivity to vari- 
ations of difl'erent parameters. The results are summarized in 
TablesElandS 

Computing the critical dynamo numbers (the threshold value 
of Cs for which the magnetic energy begins to grow), starting 
from a dipolar configuration and then from a quadrupolar one, 
enables us to test the influence of the MC amplitude and of the 
difl'usivity on the parity selection in our various cases. In a rel- 
atively low range of MC amplitude (vq < 1000 cm. s"^) and at a 
magnetic difl'usivity of 5.10^^ cm^.s"^ dipolar solutions are eas- 
ily excited in the unicellular case. On the contrary, for the mul- 
ticellular models, the symmetric parity is already appearing at 
lower values of vq. Indeed, the magnetic field has switched to 
a quadrupolar parity at vq = 785 cm.s"^ for the two radial cell 
model and the diff'erence between Cf (A) and Cf(S ) is already 
very small at vq = 643 cm.s"^ (0.84 compared to 0.85). The 4- 
cell model always favours the quadrupolar configuration, even 
for low velocities. 

To obtain a smaller Cf^^ for the dipolar mode in the multiple 
cell models, we need to increase the magnetic diff'usivity. Indeed, 
for all cases, increasing the diff'usivity widens the range for vq in 
which we stay in the dipolar configuration. In case 1, at 7/t = 
5.10^0 cm^.s-i andvo = lOOOcm.s'^ = 1 .4 for the dipole 
and Cf^^ = 1.35 for the quadrupole, which explains the drift of 
parity we observe in this case in Table [2l When we increase the 
magnetic diff'usivity up to rjt = 8.10^^cm^.s"\ the dipole be- 
comes the most easily excited solution with Cf^^(dipole) = 1.29 
and Cf^^ (quadrupole) = 1.3. In the same way, the systematic 
parity switching in the 4-cell model disappear when we go from 
T]t = 5.10^^cm^.s"^ to r]t = 8.10^^cm^.s"^ and for example at 
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Table 2. Critical values of Cs starting from a dipole (A) 
or a quadrupole (S) for various values of vq and at r/t = 
5.10^^ cm^.s"^ for the 3 configurations of meridional circulation. 
The favoured symmetry (the smallest values of Cs) is indicated 
in bold characters. In the first line the Cf^^ for the reference case, 
case la and case 2a are shown. 





Single Cell 


Two radial cells 


4 cells 


vo 








643 


1.88 1.92 


0.84 0.85 


1.16 1.01 


785 


2.34 2.37 


1.05 1.04 


1.19 1.05 


1000 


3.07 3.05 


1.40 1.35 


1.55 1.38 


1500 


4.51 4.49 


2.50 2.20 


3.00 2.73 



Vo = 785 cm.s"\ the solar symmetry is favoured. We thus con- 
firm the work of Chatterjee et al. (12QQ4b which shows that in- 
creasing the diff'usivity in the convection zone, thus allowing 
diff'usive coupling of the poloidal field between the two hemi- 
spheres improves the parity conservation. 

Table 3. Critical values of Cs starting from a dipole (A) 
or a quadrupole (S) for various values of vq and at t]i = 
8.10^^ cm^.s"^ for the 3 configurations of meridional circulation. 
The favoured symmetry (the smallest values of Cs) is indicated 
in bold characters. 





Single Cell 


Two radial cells 


4 cells 


Vo 


C^^(A) C^'(S) 


C^^(A) C^\S) 


C^^(A) C^'(S) 


643 


2.72 2.81 


0.79 0.85 


2.10 2.14 


785 


2.79 2.95 


0.99 1.03 


2.11 2.15 


1000 


2.89 3.08 


1.29 1.30 


1.80 1.50 


1500 


4.41 4.43 


1.96 1.94 


2.12 1.92 



However, as soon as we increase the amplitude of the ve- 
locity field via an increase of the MC amplitude, we recover 
the parity drift from a dipolar to a quadrupolar configuration, 
the quadrupole becomes the easiest solution to excite. Figure [TOl 
shows a typical representation of a parity shift in the case of 4 
cells of meridional circulation per hemisphere: the toroidal field 
at the base of the CZ is switching from an antisymmetric con- 
figuration with respect to the equator to a symmetric one. As 
Dikpati & Oilman ( 20011) showed, it is the connection at the 
equator of the sufficiently strong poloidal fields of each hemi- 
sphere that enables the cancellation of Bf and the creation of an 
antisymmetric Bq. The shear of this antisymmetric Bq by dif- 
ferential rotation is then responsible for the creation of anti- 
symmetric toroidal field which we observe in the Sun. Thus it 
is very likely that increasing the velocity amplitude makes the 
magnetic field travel faster along the conveyor belt and prevents 
the poloidal field from staying in the same location enough time 
to connect with its counterparts in the opposite hemisphere. As 
a consequence, models with faster flows shift to quadrupolar so- 
lution since they prevent these reconnection phenomena. 

This particular property of parity selection explains why it is 
not sufliicient to act on the MC amplitude in case 1 to recover a 
satisfying model with a 22-yr cycle. Indeed, increasing the ve- 
locity causes the magnetic field to become symmetric with re- 
spect to the equator. On the contrary, we saw that increasing 
the diflTusivity tends to favour dipolar symmetry. As a conse- 
quence, cases lb and 2b, which have relatively strong values of 
the diflTusivity are able to reproduce the field antisymmetry we 
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Time (years) 

Fig. 10. Zoom on the epoch of parity drift of the toroidal field 
for the 4-cell model with uq = 2000cm.s"\ = 20cm.s"^ and 
j]t = 5.10^^ cmls-^ 



observe in the Sun, even with significant meridional flow am- 
plitudes. For case lb, the values of the Cf are Cp(A) = 2.47 
and Cp(^) = 2.48 and for case 2b, we get Cp(A) = 3.36 and 
Cp(5') = 3.39. These cases are thus able to sustain a 22-yr cycle 
period without drifting from a dipolar to a quadrupolar configu- 
ration. 



5. Conclusion and perspectives 

In this paper, we have discussed 2D BL flux transport type solar 
dynamo models with various profiles of meridional flows. 

We have first tested the influence of introducing the parame- 
ters giving a solar-like solution in the reference unicellular case 
in the muticellular models. These cases, denoted la and 2a, 
show that the presence of a multicellular circulation has a strong 
perturbing impact on the behaviour of solar dynamo models. 
Adding cells in radius (case la) leads to a complex advective 
path and thus causes the cycle period to be more than tripled 
compared to the reference case (the cycle is here of 84.6 years 
instead of 22 years in the reference model) and in this case, a 
strong poleward branch appears on the butterfly diagram, due to 
the poleward advection by the MC at the base of the CZ. The ra- 
dial field at the surface seems to show very fine and small struc- 
tures during its whole cyclic evolution. In this model, we notice 
that the cycle period is moreover strongly linked to the ampli- 
tude of the meridional flow, which indicates that the cycle will 
be significantly sensitive to the observed fluctuations in the MC 
amplitude. In the 4-cell model (case 2a), the most obvious prop- 
erty is that we seem to get two magnetic cycles, with diflTerent 
periods, both longer than in the reference case (44.7 years near 
the equator and 124 years near the poles). Moreover, the phase 
relationship between the poloidal and toroidal parts of the field 
does not match the solar observations anymore. However, unlike 
case la, the dependance on the amplitude of the MC is reduced, 
which could make the model and thus the cyclic activity more 
robust and less sensitive to temporal fluctuations observed in the 
Sun. 

The set of parameters were in these cases clearly not adapted 
to recover a 22-yr cycle period, we thus modified the appropri- 
ate parameters to get a solar-like solution. Relying on the least- 
square fits, cases lb and 2b were thus computed with higher 
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difFusivities and higher MC amplitudes and a 22-yr cycle was 
indeed recovered. It should be noted that we stay in these cases 
in realistic values for both magnetic diffusivity and MC ampli- 
tude and that our models are still all dominated by advection, 
as shown on the evolutions of the field lines in the meridional 
plane. For these cases, we note that the butterfly diagrams are 
smoothed out probably thanks to the increase of magnetic difl'u- 
sivity, the small structures visible on the radial field on cases la 
and 2a vanish, we thus seem to obtain a globally less perturbed 
behaviour for these models. In case 2b, we moreover decrease 
the intensity of the polar surface field, which is in better agree- 
ment with observations. But we still obtain properties such as 
the phase relationship between the poloidal and toroidal parts of 
the field that are not that which we observe in the Sun. Since 
parity selection is one of the major concerns in BL flux transport 
dynamo models, we focus in Sect. [Hon the parity issue in these 
models. We show that adding cells both in radius and in latitude 
seems to favour the quadrupolar parity, which we do not observe 
in the Sun. However, if the magnetic difl'usivity is sufficiently 
high, we get diff'usive coupling of the poloidal field across the 
equator and thus the dipolar parity conservation is improved. On 
the contrary, if the MC amplitude is increased, the major trend 
of all cases is to switch from a dipolar to a quadrupolar magnetic 
field configuration. We can nevertheless recover cases with com- 
plex meridional flow with a 22-yr cycle period and a favoured 
dipolar parity, staying in realistic values of both MC amplitude 
and magnetic difl'usivity, this is the case for models lb and 2b. 

As far as the meridional flow amplitude is concerned, we 
shall note that our relatively small velocities at the solar surface 
are related to the small stratification we have in our models (our 
density is proportional to 1/r) which implies a small velocity 
contrast between the surface and the bottom of the convection 
zone. Indeed, the density profile used by Dikpati & Charbonneau 
(11999b was proportional to ^jR^lr - 1 which has a strong vari- 
ation in radius especially near the surface. Thus, the most im- 
portant velocity amplitude (that at the base of the CZ, which ad- 
vects the strong toroidal field created by diff'erential rotation) can 
be identical with very diff'erent velocities at the surface in these 
two models using diff'erent density profiles. We also note that 
since there was a factor 2 between vq and max(v6i) in the work 
of Dikpati & Charbonneau (I1999K a direct comparison with our 
work (where vq = max(v6i)) is not straightforward. 

Even if cases lb and 2b seem to allow a 22-yr cycle com- 
bined with persistent solar-like dipolar parity and a smooth 
cyclic butterfly diagram, we show that introducing a complex 
MC in our models has a strong perturbing impact, for exam- 
ple on the phase relationship between the poloidal and toroidal 
parts of the magnetic field which does not correspond to the so- 
lar observations. We are thus heading to the hypothesis that the 
BL mechanism may not be the only source of poloidal field in 
the solar dynamo cycle, especially if the Sun happens to show a 
persistent multicellular meridional flow which seems to be quite 
destructive for several solar cycle features in the pure BL flux- 
transport framework. Of course we now need to check the influ- 
ence of a less monolithic meridional flow structure, with extra 
cells more concentrated in a particular area of the CZ and vary- 
ing in time since the position and strength of each meridional 
cell seem to influence quite significantly the global properties of 
the solar dynamo. 

It is thus now a real challenge for local helioseismology to 
probe the Sun deeper to give better constraints on the meridional 
flow in the convection zone. 
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Appendix A: Numerical approach 

To solve the equations, we use a code adapted by P. Charbonneau 
and T. Emonet from Finite Element Analysis by D.S.Burnett 
(119871) . This code enables us to solve a general partial differen- 
tial equation (PDE) using a finite element method in space and 
a third order scheme in time. We adapted it to problems such as 
Qf^ - Q, flux transport or multicellular flux transport solar dy- 
namos and we implemented new boundary conditions (radial or 
potential field at the top) and initial conditions. 



A.1. Spatial method 

The finite element method is a very eflScient way to obtain ap- 
proximate solutions to linear or non linear PDEs in any kind of 
geometry. Our code STELEM (STellar ELEMents) solves Eqs.[T] 
and |2] with this method, ie seeking the approximate solutions 
and B(p as linear combinations of trial functions i/^i (to be more 
specific these trial functions are Lagrange polynomials of degree 
1 (linear functions) and serendipity shape functions for second 
order interpolation (quadratic functions), depending on the com- 
plexity of our equations). 



The main steps of the method are the following: 

Our domain (the annular meridional cut) is divided 
into smaller regions called elements. In our case, they 
are straight- sided quadrilaterals when we use first order 
Lagrange polynomials, with a node at each corner of the 
quadrilateral and with a node at each corner and one extra 
node per side and without any interior node if we use second 
order interpolation (see Fig. lA.lb . 

In each element, the PDEs are transformed into ordinary dif- 
ferential equations (ODE) in time involving the coefficients 
ai{t) and biif) of the linear combinations. 
The terms in the element equations are numerically evalu- 
ated for each element in the mesh. The resulting numbers 
are assembled into a much larger set of equations called the 
system equations. 

The boundary conditions are taken into account. They can 
be of Dirichlet type (we impose the value of the function 
at the boundary) or of Neumann type (we impose the nor- 
mal derivative of the function at the boundary). In particu- 
lar for the potential extrapolation, we proceed as the follow- 
ing, we were largely inspired by the procedure of Dikpati & 
Choudhuri[l99l 

The top boundary condition is that we have to match 
smoothly our magnetic field B(r, 6, t) to a field satisfying the 
free space equation: 

V X B(r, 6/, = 

As we work in spherical axisymmetric geometry, we write 
that: 

B(r, 6/, = V X (A^(r, 6, t)e^) -h B^(r, 6>, Oe^ 

And the equation of free space leads to two equations, one 
concerning ^^(r, 6, t) and the other one concerning A0(r, 6, t), 
which are: 




Bilinear element 



Biquadratic element 



Fig. A.l. Sketch of the quadrilateral mesh we are using in the 
meridian plane, uniformly spaced in r and more accurate in the 
polar regions. As we work with the variables cos 0,r, we get a 
rectangular grid in the r,cos 6 plane. On the right, we show a 
zoom on a single rectangular cell with the 4 nodes at each corner 
in the case of the first order interpolation and with one extra node 
per side in the case of the second order interpolation. Note that 
the cells in the quadratic case do not contain any interior node. 



06 



dr 



(V^ 



1 



sin^ 6 



= 







As we are dealing with a finite element method, the most 
convenient and natural procedure is to seek to express these 
boundary conditions as either Dirichlet or Neumann condi- 
tions for A^ and 5^. 

Equation for 5^ very easily leads to the solution 5^ = 
C/(r sin 6), C being a constant. We fix this constant value to 
so that our top boundary condition on 5^ is the homogenous 
Dirichlet condition: B^{Rq, 0, t) = 0. 

We now have to deal with the more difficult condition on A^. 
A general solution to this equation can be written in the form 



A^(r>RQ,0,t) = Y,^Pi(cose) 



where P^icosO) is the associated Legendre polynomial. We 
find that truncating the sum at Ne/2, half of the number of 
grid points in 6 always result in an error in the projection of 
less than 10"^. Hence, the coeflftcients an(t) are the coeffi- 
cients of the expansion of 7?"^^A0(7?q, 6, t) on the associated 
Legendre polynomials. Thus, the value of Unit) is calculated 
by the scalar product of R^^^A^{Rq, 6, t) with P^(cos^), di- 
vided by the norm of leading to: 



an{t) 



^ M^e^ ^' O^i(cos 0) sin OdO 
£ [Plicos 6)]^ sin OdO 



By the variable change x = cos(6) in the upper integral, we 
are led to calculate the integral on [-1,1] of the product of 



Jouve Laurene & Brun Allan Sacha: On the role of meridional flows in flux transport dynamo models 



13 



two smooth functions. Thus we calculate this integral using a 
Gauss-Chebyshev quadrature formula which uses the weigh- 
ing functions 1 / Vl - x^. The lower integral is the norm of 
the associated Legendre polynomials, we know that its value 
is: 2n{n+ 1)1 (2n+ 1). 

From the coefficients a„(0, we can deduce the derivative of 
at the solar surface: 

n-l O 

and from a simple finite difference scheme, we impose a 
Dirichlet condition on the poloidal potential, calculating the 
new value of at the surface, using the points of the layer 
immediately below the surface. It leads to: 



A^(/?0, e, t) = A^Rq - Ar, e, t) + Ar— 1,.^^ (A.l) 
= A,(Re - Ar, 0, - Ar ^ ^Jl^^^P^^cos 6) 

Inside the same time step, we have then to recalculate the co- 
efficients Unit) with the new value of A0(7?q, 6, t) until we get 
a sufficiently small difference between two successive values 
of A0(7?Q, 6, t) in order to make the procedure converge, we 
usually do not need more than 10 iterations to get a conver- 
gence with a relative error of 10"^. 
- We get a final set of ODEs in time which we solve with a 
third order scheme we describe below. 



A.2. Temporal scheme 

The scheme that we use is adapted from Spalart et al. 119911 We 
have to solve the following ODE: 

N being the non linear operator evaluated in the preceeding 
step (the finite element method). 

The three steps of this explicit scheme enable us to get an 
error as small as o(At^), the different steps are: 

If Un = u(t) and = u(t -\- At), it leads to: 

u = Un-\- AtyiN(un) 

u' = u + At[y2N(u) -h ^iN(un)] 

Un+l = -h At[y3N(u') -h ^2N(u)] 

The coefficients yi, 72, 73, ^1 et ^2 are deduced from the 
Taylor expansion of u(t -h AO and thus leads to: 71 = 8/15, 72 = 
5/12, 73 = 3/4, = -17/60 and ^2 = -5/12. 



A.3. Code validation 

The STELEM code was validated thanks to an international dy- 
namo benchmark in Jouve et al. 120071 

All data and notes can be fo und at the address: 
http://www.nordita.dk~brandenb/tmp/benchmark 



